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Abstract 

Bombus flavescens Smith is one of the most widespread bumblebee species in the Oriental region. Due to 
colour polymorphisms, this species or species-complex has been a challenge for taxonomy. This study aims 
to assess the taxonomic status of the flavescens-complex using evidence from COI barcodes and morphol- 
ogy. We then reconstruct its biogeographic history from a phylogenetic analysis of populations across the 
current range, combining COI with 16S and nuclear PEPCK data. Despite a large range of polymor- 
phisms across its distribution, the results show that B. flavescens is a single species based on algorithmic 
species delimitation methods, and it is clearly separated from its sister species, B. rotundiceps Friese. We 
suggest that B. flavescens diverged from its sister lineage in the Himalaya and dispersed into Southeast Asia 
in the Pleistocene. Conservation of the widespread B. flavescens will need to consider its several unique 


island populations. 
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Introduction 


Bumblebees (genus Bombus Latreille) are well-studied pollinators, especially in tem- 
perate regions (Williams et al. 2014; Rasmont et al. 2021). They can also be found in 
both subtropical and tropical areas in Central and South America, and Asia (Williams 
1998), although taxonomic knowledge gaps and the relative rarity of bumblebees in 
these areas still constrain our understanding of their biogeography and ecology. 

The taxonomic studies of bumblebees in Asia, based on morphological evidence 
only, have been considered highly problematic, due to colour pattern polymorphism 
within the same species (Huang et al. 2015b; Ding et al. 2019) and cryptic species 
(Williams et al. 2020, 2022a). Moreover, the same colour pattern is also observed from 
different species which are co-existed in same locality (Hines and Williams 2012). 
Therefore, these confusions of taxonomic delimitation of Asian bumblebees require 
molecular evidence (Williams et al. 2012). Taxonomic status of bumblebees has been 
revised and various new species have been recently described in Asia, based both mor- 
phology and molecular evidence (Williams et al. 2020, 2022a, c; Williams 2022). 

Bombus flavescens Smith is a Southern Asian species of the subgenus Pyrobombus 
Dalla Torre (Williams 1998). While the type locality of B. flavescens is Zhoushan, 
China (Smith 1852), this taxon is found in many other countries in the Himalaya and 
Southeast Asia with tropical low montane habitats, including Nepal, India, Bhutan, 
Myanmar, Thailand, Laos, Vietnam, Malaysia, and the Philippines (Starr 1992; Wil- 
liams 1998, 2022; Williams et al. 2009; Koch and General 2019). 

Most frequently, B. flavescens is a predominantly black pubescence (or hair) bum- 
blebee with an orange tail (the area of hairs covering posterior part of metasoma) and 
legs in workers and queens. Males often (e.g., in parts of China and in the Himalaya) 
show a predominantly pale yellow (flavescent) hair pattern with orange legs (Fig. 1), 
although they sometimes exhibit the widespread dark pattern of the females. However, 
the colour pattern of hair varies across its wide area of distribution (Fig. 2), which has 
caused taxonomic problems for more than a century and led to the formal descrip- 
tion of local variants (Frison 1934). For example, the name B. alienus Smith in China 
(Fig. 2E; morphologically synonymised by Williams (2022)) was applied to a specimen 
with an anterior yellow band on its metasoma (terga 1—2 or T1—2) whereas specimens 
of from the Philippines (Fig. 2G), with yellow hair covering the side of the mesosoma 
and metasoma (T1—2) but without the orange tail, have been described as B. bagui- 
onensis Cockerell (morphologically synonymised by Williams (1998)). In addition, a 
form with an entirely body, covered with uniformly orange hairs, has been described as 
B. rufoflavus Pendlebury (Fig. 2J; provisionally synonymised by Williams (1998), based 
on morphology). Further names were published for colour patterns resembling the typi- 
cal B. flavescens: B. mearnsi Ashmead from Mindanao in the Philippines (morphologi- 
cally synonymised by Pittioni (1949)); B. bakeri Cockerell from Negros (Fig. 2H; mor- 
phologically synonymised by Frison (1925)); and B. tahanensis Pendlebury from the 
Tahan mountains of Malaysia (Fig. 21; morphologically synonymised by Frison (1934)). 
Apart from this high level of colour variation, other morphological characters, includ- 
ing shape of labrum, punctures on clypeus, and male genitalia, are relatively similar 
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Figure |. Holotype male of B. flavescens, deposited in the entomological type collection at the Natural 


History Museum, London (NHMUK 014025381). Scale bar: 5 mm. 
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Figure 2. Colour variation of B. flavescens workers with their unique specimen identifiers A India (FLA#5) 


B Nepal (NHMUK 010814958) C Thailand (CT#1024) D Kunming, China (FLA#6) E alienus; Fujian, 
China (NHMUK 010815023) F Taiwan (FLA#4) G baguionensis; Luzon, the Philippines (CT#932) 
H bakeri; Negros, the Philippines (ZMA.INS.758600) I tahanensis, Gunung Tahan, Malaysia (NHMUK 
010820651) J rufoflavus; Cameron highlands, Malaysia (NHMUK 010814574). Some images are re- 


versed. Scale bar: 5 mm. 


(Williams 1998; Thanoosing 2022). Yet, at various times through the 20" century dif- 
ferent intraspecific names have been applied to describe minor B. flavescens variants (Fri- 
son 1925, 1928, 1934; Bischoff 1936; Chiu 1948; Pittioni 1949; Table 1) although they 
are clearly grouped together because of their morphological similarity (Williams 1998). 
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Table |. List of Bombus flavescens names. 


Author Name 
Smith (1852) Bombus flavescens 
Smith (1854) Bombus alienus 
Ashmead (1905) Bombus mearnsi 
Friese (1905) Bombus rufocaudatus 
Cockerell (1917) Bombus geei 
Cockerell (1920) Bombus irisanensis vat. baguionensis, Bombus bakeri 
Pendlebury (1923) | Bombus tahanensis, Bombus rufoflavus 
Frison (1925) Bremus mearnsi, Bremus mearnsi var. bakeri, Bremus irisanensis vat. baguionensis 
Hedicke (1926) Bombus imuganensis 
Frison (1928) Bremus (Pratobombus) baguionensis, Bremus (Pratobombus) baguionensis vat. imuganensis 
Dover (1929) Bremus rufoflavus 
Skorikov (1933) Pratobombus flavescens 
Frison (1934) Bremus mearnsi vat. deflectus, Bremus mearnsi vat. ditutus, Bremus mearnsi vat. bakeri, 


Bremus mearnsi vat. geei 

Bischoff (1936) Bombus (Pratobombus) mearnsi ssp. chekiangensis 

Chiu (1948) Bombus mearnsi vat. deflectus, Bombus mearnsi vat. ditutus, Bombus mearnsi vat. bakeri, Bombus 
mearnsi vat. geei, Bombus mearnsi var. subrufus, Bombus mearnsi var. luteus 

Pittioni (1949) Bombus flavescens f. dilutior 


Colour-pattern variation of bumblebees can be controlled by differential gene ex- 
pression of, e.g., the Abd-B and nubbin genes affecting pigment generation during pu- 
pal development (Rahman et al. 2021). However, the evolutionary forces determining 
variation in colour pattern among B. flavescens populations remain unclear. Selection 
for mimicry might take a major role in driving evolutionary differences because the 
colour pattern of B. flavescens locally matches other bumblebee species in the same 
areas (Williams 2007). For example, on Luzon Island, baguionensis shows a similar pat- 
tern to B. (Megabombus) irisanensis Cockerell, and in the Malay Peninsula, rufoflavus 
resembles B. (Megabombus) montivagus Smith with a similar predominantly orange 
pattern of both pubescence and underlying sclerites, which had been named B. max- 
welli Pendlebury (Hines and Williams 2012). 

Orange body colour (pubescence and sclerites) and enlarged ocelli in Hymenoptera 
are often associated with a nocturnal or crepuscular lifestyle, for example, Megalopta 
bees and Apoica wasps in Central and South America (Roubik 1989; Williams 2007; 
Warrant 2008). Interestingly, male and female nocturnal carpenter bees, Xylocopa 
(Nyctomelitta) myops Ritsema, show a pattern closely similar to the orange pattern 
bumblebees. ‘This carpenter bee is found in Malaysia, Singapore, and Borneo, based on 
specimens in the Natural History Museum, London (NHMUK) collection and records 
by Ascher et al. (2022). Similarly, the two orange pattern bumblebees in Malaysia, 
rufoflavus (B. flavescens) and maxwelli (B. montivagus), might be active nocturnally 
also, at least in part (Williams 2007). Species of the subgenus Nyctomelitta Cockerell 
display relatively large ocelli compared to day-flying carpenter bees (Cockerell 1929; 
Michener 2007). However, the ocelli of the orange pattern bumblebees remain to be 
assessed in this regard. Observation of a night-flying bumblebee had been claimed 
in Myanmar (Doria 1886), but this record is a misidentification of the carpenter 
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bee, X. (Nyctomelitta) tranquebarica (Fabricius) (Cameron 1910). Therefore, genuine 
evidence of nocturnal bumblebees has yet to be recorded in Southeast Asia. 

The colour-pattern diversity of B. flavescens must be seen in light of the complex bioge- 
ography of its distribution range across the Southeast Asian mainland and the Philippines. 
The most recent common ancestor (MRCA) of the flavescens species complex has been 
placed into the Miocene epoch around eight million years ago (Ma), presumably on the 
mainland (Hines 2008) with two possible dispersal routes to the Philippines (Starr 1989): 
1) through the Sundaland route, via Borneo—Sulawesi from the South; or 2) across the 
Luzon strait via Taiwan from the North. When global temperatures increased after the last 
glacial maximum (LGM), bumblebees likely became restricted to highlands (e.g., the Cam- 
eron Highlands, Malaysia) and diverged in allopatry. At the same time, as the sea level rose 
and submerged land bridges across Sundaland (Voris 2000; Sathiamurthy and Voris 2006; 
Woodruff 2010), populations of bumblebees on the Philippines were even further isolated. 

This study seeks to clarify genetic patterns within the flavescens-complex, to estab- 
lish the species status of geographically separated lineages and to infer biogeographic 
scenarios for bumblebees in Southeast Asia more generally. Species-level entities in 
bumblebees have been delimited by morphological, molecular (e.g., Williams et al. 
(2019, 2020)), and chemical approaches, the latter using Cephalic Labial Gland Se- 
cretions or CLGS (e.g., Brasero et al. (2021), Ghisbain et al. (2021)). Given the chal- 
lenges of morphology for taxonomic assessment of the flavescens-complex, molecular 
data are required to resolve the species status of sub-lineages and their relationships. In 
addition, this study aims to clarify the question about the nocturnal lifestyle in B. fla- 
vescens, using morphological examination of their ocelli. 


Methods 
Sampling 


Museum and institutional specimens of the flavescens-complex (Table 1) were exam- 
ined for morphological analysis. Collecting information was recorded and made avail- 
able at the NHMUK Data Portal (https://doi.org/10.5519/qd7f4uuw and https://doi. 
org/10.5519/isxh6saw). Additional records were obtained from reliable sources, in- 
cluding: 1) published literature (Williams et al. 2009, 2010); 2) major natural history 
collections, including the Naturalis Biodiversity Center (NMNL; Bakker F, Creuwels 
J), the Smithsonian Institution (USNM; Orrell T, Informatics Office), the University 
of Illinois at Urbana-Champaign (INHS; McElrath T), and the Taiwan Forestry Re- 
search Institute (TFRI; Lu S), available on the Global Biodiversity Information Fa- 
cility (GBIF) (GBIEorg 2021). Records with unclear locality information and clear 
geographic outliers were ignored. For example, there is a record of B. flavescens from 
Japan (e.g., GBIF occurrence 3801818589; McElrath 2022), far outside the estab- 
lished range of B. flavescens, based on the specimen records by Williams (1998, 2022). 

COI barcode data of B. flavescens and relatives were obtained from public databases 
(Suppl. material 1) and newly sequenced from pinned or fresh specimens (Table 2). 
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Table 2. List of specimens included, species name, ID, deposit place (KKIC = Kasetsart Kamphaeng Sean 
Insect Collection, PHW = Paul H. Williams Research Collection, CUNHM = Chulalongkorn Univer- 
sity Natural History Museum Collection, NMNL = Naturalis Biodiversity Center Collection), sex/caste 


(w = worker, m = male), locality, collecting date, and COI GenBank accession number. 


Species ProjectID Collection (specimen ID) Sex/caste Locality Collecting GenBank 
date accession 
B. flavescens CT#662 KKIC w Thailand, Nakorn = 7/5/2017 OP355718 
Pathom? 
CT#669 KKIC m Thailand, Loei 13/4/2016 OP355719 
CT#926 PHW w Malaysia, NA OP355720 
Gunung Tahan 
CT#1018 NMNL (ZMA.INS.758598) w Philippines, Negros 4/5/1953. OP355722 
CT#1023 CUNHM (BSRU-AB-9609) w Thailand, 15/2/2021 .OP355725 
Chiang Mai 
CT#1024 CUNHM (BSRU-AB-9610) Ww Thailand, 15/2/2021 OP355724 
Chiang Mai 
FLA#2 PHW w Philippines, Luzon 27/4/1986 OP355725 
FLA#3 PHW Ww Taiwan, Tai Chung 5/5/1980 OP355726 
FLA#4 PHW w Taiwan, Nantou 24/6/1989 OP355727 
FLA#5 PHW w India, Uttar Pradesh 6/5/1990 OP355728 
FLA#6 PHW w China, Kunming 8/4/2018 OP355729 
FLA#7. NMNL (ZMA.INS.773029) Ww Bhutan, Lungtenphu 10/6/1996 OP355730 
FLA#10 PHW Ww Philippines, Luzon 27/4/1986 OP355731 
B. rotundiceps CT#960 CUNHM (BSRU-AB-1268) m Thailand, 17/6/2019 OP355721 
Chiang Mai 
ROT#1 PHW w India, Uttar Pradash 6/5/1990 OP355732 
ROT#2 PHW Ww India, Uttar Pradash 6/5/1990 OP355733 
ROT#3 PHW Ww China, Guangxi 3/6/2016 OP355734 
ROT#4 PHW w China, Guangxi 5/6/2016 OP355735 


The samples were selected to represent the geographical distribution and all major col- 
our pattern variants of B. flavescens (Fig. 3). Bombus flavescens belongs to the pratorum- 
group (Williams 1998). Thus, closely related species in the pratorum-group (Cam- 
eron et al. 2007; Williams et al. 2009, 2010), were included in the dataset: B. ardens 
Smith, B. pratorum (Linnaeus), B. pyrenaeus Pérez, B. modestus Eversmann, B. nursei 
Friese, B. biroi Vogt, B. wangae Williams et al., and B. rotundiceps Friese. The status of 
B. nursei as a distinct species was suggested in Williams (2022). Thirty COI sequences 
(> 600 bp) of B. flavescens and the relatives were downloaded from the Barcode of Life 
Data System (BOLD; Ratnasingham and Hebert 2007) and the National Center for 
Biotechnology Information (NCBI) database (GenBank; Sayers et al. 2022) (Suppl. 
material 1). Nine B. flavescens sequences from China were provided by co-author MO, 
sequenced from freshly collected specimens. Thirteen further B. flavescens and five 
B. rotundiceps specimens from museum and institution collections were selected for 
sequencing (Table 2). The identifier numbers were applied to the studied specimens: 
1) “CT #n” for Southeast Asian bumblebee specimens 2) “FLA#n” and “ROT#n” refer 
to B. flavescens and B. rotundiceps specimens from outside Southeast Asia, respectively. 
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Figure 3. Map showing the distribution of B. flavescens with the pattern of hair colour in the dorsal view. 
The blue dots represent the records from museum collections, literature, and GBIF database (n = 702). 


The red dots represent the localities of barcoded specimens. 


DNA extraction 


Pinned and fresh specimens of B. flavescens and B. rotundiceps were selected for DNA 
extraction. For the recent specimens (< 20 years old), the tissue sources were a right 
front leg or a whole body for high DNA concentration yield. Each leg sample was 
ground using a small pestle in a microcentrifuge tube, whereas the whole-body sample 
was directly put in a microcentrifuge tube without specimen damages. ‘The tissue sam- 
ples were incubated at 56 °C for 24 hours in the ATL buffer with Proteinase K enzyme. 
Genomic DNA was extracted using the Qiagen DNeasy Blood and Tissue Kit, fol- 
lowing the kit protocol. DNA quality and quantity were assessed using the Nanodrop 
spectrophotometer (Thermo Scientific ND8000) and the Agilent 2200 TapeStation 
system (Agilent Technologies, Inc.). 

For older specimens (> = 20 years old), genomic DNA was likely to be degraded 
and all DNA extraction steps were performed inside a laminar flow chamber using 
dedicated UV sterilised equipment and consumables different from those used with 
modern bumblebee DNA. The DNA extraction method followed Sproul and Maddi- 
son (2017) using the Qiagen QIAamp Micro Kit, but without the use of carrier RNA. 
For each specimen, the right front leg was cleaned with sterile water, frozen with liquid 
nitrogen, and then ground with a small pestle. The samples were incubated in ATL 
buffer with Proteinase K at 55 °C for 24 hours. A blank sample (only the ATL buffer 
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with the enzyme) was included for each extraction batch as a contamination control. 
Samples were eluted from the column twice with 15 yl AE buffer and incubated at 
room temperature for 10 minutes. The DNA concentration was measured using the 
Qubit fluorometer high sensitivity (Invitrogen) and samples were stored at -20 °C. 


Primer design, amplification, and sequencing 


For the recent specimens, PCR was performed to amplify the full length of COI ampli- 
cons (658 bp) using primers LepF1 and LepR1 (Hebert et al. 2004) and an annealing 
temperature at 50 °C. For older specimens, due to the DNA degradation, three pairs of 
primers were newly designed for amplifying partial COI contigs (150-300 bp), using 
Primer3 software (Untergasser et al. 2012). The COI sequence of B. flavescens (Gen- 
Bank accession number: GU085209) was used as a reference sequence for the first two 
pairs of primers. ‘The third pair was designed based on the partial mitogenome of B. pra- 
torum (GenBank accession number: K 1164684). There are three pairs of new primers 
for partial COI amplification (contig) in this study (Table 3): 1) FLA1: ParCOIl_246_ 
FLA“ Io Cim? “= -66,2°°C) and, ParCOl “416 FLA Rh (im? =°63,7.°C), 2)-FLA2: 
ParCOI_349 FLA F2 (Tm® = 64.2 °C) and ParCOI_662_FLA_R2 (Tm° = 66.6 °C), 
and 3) PRAI: ParCOI_1816_PRA_F1 (Im° = 61.7 °C) and ParCOI_1994_PRA_ RI 
(TIm® = 65.5 °C). These primers were tested for amplification efficiency (Suppl. materi- 
als 2-4). The annealing temperature for the pair FLA1 and FLA2 was set at 65 °C and 
PRA was set at 63 °C. 

The PCR products were sequenced in both forward and reverse directions using 
ABI technology at the NHMUK’s sequencing facility. The sequences were edited using 
MEGA version 7.0.26 (Kumar et al. 2016: https://www.megasoftware.net/), to trim 
low-quality bases near the ends of the traces. The COI barcodes have been deposited 
in GenBank (Table 2). 


Phylogenetic analysis and species delimitation 


The aligned dataset was analysed to determine the best-fitting nucleotide substitution 
model using jModelTest software version 2.1.6 (Darriba et al. 2012). Phylogenetic 
trees were constructed using Bayesian Inference (BI) with MrBayes version 3.2.2 (Ron- 


Table 3. Newly designed primers were used in this study, including gene, primer name, strand (F = for- 


ward, R = reverse) and primer sequence. 


Gene Primer name Strand Primer sequence 
Partial COI ParCOI_246_FLA_F1 lf 5'— CCT GACATAGCPITEGCCACGA -3' 
ParCOI_416_FLA_R1 5’—- TGCAATATCAACTGAAGGTGATG -3’ 
ParCOI_349_FLA_F2 5'- CAGGATGAACTGTT FACCETCCT =3° 
ParCOI_662_FLA_R2 5’ TGGATCACCT CCT CCTATTGGA.=3' 
ParCOI_1816_PRA_F1 5’- TTCGTATAGAATTAAGTCATCCTGGT -3’ 
ParCOI_1994_PRA_R1 5’- TCGTGGAAAAGCTATATCAGGTGAT -3’ 


A mA TA 
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quist and Huelsenback 2003). Markov Chain Monte Carlo (MCMC) chains were 
run for 10 million generations and sampled every 1000" generation. The first 10% 
was discarded as burn-in. The phylogenetic trees were visualised using FigTree ver- 
sion 1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/). Bumblebee species from closely 
related subgenera (Cameron et al. 2007) were used as outgroups, including B. (Bom- 
bus) terrestris (Linnaeus) and B. (Alpinobombus) polaris Curtis. The outgroup sequences 
were generated by Thanoosing (2017). 

Poisson Tree Process (PTP) analysis is used to identify likely species coalescents. 
The test establishes the transition from long branches defining between-species diversi- 
fication to short branches within species (estimated from the number of substitutions 
on branches) (Zhang et al. 2013; Kapli et al. 2017). The proportion of between- and 
within-species sampling was in a range where the PTP analysis performs well (Luo et 
al. 2018). The PTP analysis was run on the bPTP server (https://species.h-its.org; ac- 
cessed 2021) under default parameters. The analysis was restricted to unique haplotype 
sequences in order to avoid introducing misleading ‘zero-length’ branches. ‘The in- 
put phylogenetic tree was estimated using MrBayes rooted with Bombus terrestris. Re- 
sults from the PTP analysis were compared to the Generalised Mixed Yule Coalescent 
(GMYC), which uses the branching rate rather than the number of nucleotide changes 
to establish the species-to-population transition (Pons et al. 2006). The required ultra- 
metric tree was obtained with BEAST (Bouckaert et al. 2019) using the same dataset 
as for the PTP analysis. The clock rate was set at 0.0177 using the divergence rate from 
tenebrionid beetle COI (Papadopoulou et al. 2010). The GMYC analysis was run in 
the splits package of R (Ezard et al. 2009; https://r-forge.r-project.org/projects/splits/). 


Biogeographic analysis 


To reconstruct phylogenetic trees for biogeographic scenarios of B. flavescens, more 
genetic markers were added to the dataset, including mitochondrial 16S rDNA (16S) 
and nuclear phosphoenolpyruvate carboxykinase (PEPCK). Due to a lack of fresh 
specimens, most of 16S and PEPCK sequences for relatives of B. flavescens were re- 
trieved from GenBank, especially from the dataset compiled by Cameron et al. (2007) 
(Suppl. material 5). 16S and PEPCK were not available from existing sources for three 
species; B. rotundiceps, B. wangae, and B. nursei. Bombus nursei is rare in inaccessible 
areas of the western Himalaya (Williams, 1991, 2022), and B. wangae also is a rare 
bumblebee from Sichuan (Williams et al. 2009). For B. rotundiceps and additional 
B. flavescens), the 16S and PEPCK were amplified using primers 16SWb (Dowton and 
Austin 1994) /874—16SIR (Cameron et al. 1992), and FHv4/RHv4 (Cameron et al. 
2007) and sequenced in both directions (Table 4). 

An ultrametric tree was constructed based on the combined dataset of COI 
(unique haplotypes), 16S and PEPCK (Suppl. material 6), using *BEAST (Ogilvie et 
al. 2017). If not available from the same individual, sequences in these three datasets 
were concatenated for different representative from the same area of distribution. 
For PEPCK, the sequence was partitioned for introns and exons. The closely related 
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Table 4. List of 16S and PEPCK sequences, included species, ID, and the 16S and PEPCK GenBank 


accession numbers, generated in this study. 


Taxa Project ID 16S PEPCK 
B. flavescens CT#926 OP354521 - 
CT#1023 OP354523 OP382631 
FLA#6 OP354524 OP382632 
B. rotundiceps CT#960 OP354522 OP382630 


B. (Pyrobombus) lepidus Skorikov, based on Cameron et al.’s (2007) tree, was used 
as an outgroup. The nucleotide substitution model for each fragment was selected 
according to the BIC criterion. The BEAUti software (Bouckaert et al. 2019) was used 
to generate an input file for BEAST. Trees were estimated under a strict clock and a 
Calibrated Yule Model as a species tree prior with the birth rate prior drawn from a 
Gamma distribution. Due to a lack of fossil evidence for the subgenus Pyrobombus, 
the clock was calibrated from molecular rate estimates based on five genes by Hines 
(2008), placing the origin to approximately 8.5 +/- 1 Ma. Accordingly, the prior 
distribution was 10.1—6.86 Ma (5% tails). The MCMC algorithm was run for 500 
million generations and sampled every 5000" generation. The trace file was visualilsed 
using Tracer (Rambaut et al. 2018). A 20% burn-in was discarded using TreeAnnotator 
(Bouckaert et al. 2019) to construct the maximum clade credibility tree. 

S-DIVA analysis (Yu et al. 2010) and BioGEOBEARS with DIVALIKE+] model 
(Matzke 2013a, 2013b, 2014) in the RASP package (Yu et al. 2020) were used to 
estimate biogeographic scenarios for the flavescens-complex and their relatives. The bio- 
geographical distributions were divided into 11 principal areas, based on bumblebee 
global biogeography (Williams 1996), areas of endemism in Southeast Asia and the 
distributions of flavescens-complex and their relatives (Table 5): 1) Europe (A), 2) Cen- 
tral Asia (B), 3) North and Central China and North Asia (C), 4) Korean peninsula 
and Japan (D), 5) South China (E), 6) Himalaya (F), 7) Thailand (G), 8) Peninsular 
Malaysia (H), 9) Taiwan (I), 10) Luzon (J), and 11) Negros (K). These areas were cho- 
sen so that the maximum range size for any species was limited to three of these areas, 
due to acknowledged problems with the methods in inheriting broader distributions 
(Lamm and Redelings 2009). ‘Then, the dispersal-corridor model was created with pos- 
sible bumblebee permitted dispersal routes (Fig. 4). According to this model, the pos- 
sible dispersal areas in this study are AB, ABC, ABE, AC, ACD, ACE, BEF, BE BFG, 
CD, CDE, CE, CEE CEG, CEI, EK EFG, EFI, EG, EGH, EGI, EI, EJ, FG, FGH, 
GH, GHK, HJK, HK, JJ, IJK, and JK. The burnin was set at 20%. Then, the *BEAST 


trees were used as the input trees. 


Morphological study of ocelli in the flavescens-complex 


The idea of the nocturnal lifestyle in the orange pattern bumblebees, rufoflavus (B. fla- 


vescens) and maxwelli (B. montivagus) from Peninsular Malaysia, has been introduced by 


The oriental bumblebee Bombus flavescens Say 


Table 5. Principal areas of endemic distribution of the flavescens-complex and their relatives in this study. 


Area Principal ranges included Species recorded 
A Europe B. pratorum, B. pyrenaeus 
B Central Asia B. biroi 
C North China, and North Asia B. modestus 
D Korea peninsula and Japan B. ardens, B. modestus 
ke South China B. flavescens, B. lepidus, B. rotundiceps, B. wangae 
F Himalaya B. flavescens, B. lepidus, B. nursei, B. rotundiceps 
G Thailand B. flavescens, B. rotundiceps 
H Peninsular Malaysia rufoflavus (B. flavescens), tahanensis (B. flavescens) 
I Taiwan B. flavescens 
J Luzon baguionensis (B. flavescens) 
K Negros bakeri (B. flavescens) 


C 
North Asia/ 


Central Asia 


Figure 4. A corridor-dispersal model diagram of the flavescens-complex and their relatives in this study. 


Williams (2007). Although orange hair pattern can be recognised in diurnal bumblebees 
(e.g., B. humillis liger, B. pascuorum (Scopoli), and B. muscorum (Linnaeus)), rufoflavus 
and maxwelli show the uniformly orange pattern both of pubescence and sclerites unique- 
ly. This orange pattern can be recognised in the female specimens of nocturnal carpenter 
bees, X. myops, mentioned in Williams (2007), which can be also found in Peninsular Ma- 
laysia. For this reason, we selected the female specimens of X. myops and their related spe- 
cies, X. tranquebarica from the NHMUK collection, included in this study to test the idea. 

Due to the sexual dimorphism bias and the availability of male specimens, only fe- 
male or worker specimens were chosen in this study. Female flavescens-complex specimens 
(n = 107) together with its relative, B. rotundiceps (n = 12), outgroups including the diurnal 
B. irisanensis (n = 4), the orange pattern B. montivagus (taxon maxwelli) (n = 3), and the 
nocturnal carpenter bees, X. myops (n = 20) and_X. tranquebarica (n = 12), were selected for 
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Figure 5. The measurement of morphological characters A dorsal aspect: intertegular distance (ID) 
B frontal aspect: median ocellus width (MOW) C frontal aspect: head width (HW). The specimen shown 
is B. flavescens, NHMUK 010814574, from the Cameron Highlands, Malaysia. 


morphological character measurement, including median ocellus width (MOW), inter- 
tegular distance (ID), and head width (HW) (Fig. 5), under a Wild Heerbrugg M4A ster- 
eomicroscope, calibrated with ocular and stage micrometres. The specimens of the flaves- 
cens-complex were divided into eight groups: 1) Thailand (n = 4), 2) baguionensis (n = 15), 
3) Himalaya (n = 16), 4) China (n = 20), 5) Taiwan (n = 6), 6) bakeri (n = 6), 7) rufoflavus 
(n = 20), and 8) tahanensis (n = 20). The data are available in the Suppl. material 8. 

The morphological values, including MOW, the ratio between MOW and ID 
(MOW-:ID), and the ratio between MOW and HW (MOW:HW), were visualized in 
box plots using R package ggplot2 (Wichham 2016). To test for differences between 
species, and within the flavescens-complex, statistical tests were performed in R (R core 
team 2021). First, a Shapiro-Wilk normality test was used to examine the distribution 
of data of each of the flavescens-complex taxa and each bee species. Where the data 
showed a normal distribution, an ANOVA test and Tukey test were then performed. A 
Kruskal-Wallis test and Dunn’s test were used instead if the data were not normally dis- 
tributed. The Dunn’s test was conducted using the R package FSA (Ogle et al. 2022). 


Results 


DNA extraction, amplification, and sequencing 


Genomic DNA was extracted successfully for most specimens of the flavescens-complex 
and B. rotundiceps. Most of the old samples were degraded and had relatively low DNA 
concentrations (< 10 ng/ul). The three amplicons obtained with new primer pairs PRA1, 
FLA1, and FLA2 were assembled into a single contig. However, there was a 25-base pairs 
gap between the first (PRA1) and second (FLA1) fragments, and the 5’ part of the COI 
barcode region (88 base pairs) was not amplified at all. COI sequences were generated 
for 18 collected specimens with top BLAST hits to the subgenus Pyrobombus. 
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Phylogenetic analysis and species delimitation 


The final dataset of the flavescens-complex and relatives comprised 57 COI sequences, 
including 25 unique haplotypes. ‘The resulting 575 bp alignment was used to construct 
the BI phylogenic tree under the GTR+I° model selected by jModelTest (BIC = 4596) 
(Fig. 6). Bombus rotundiceps, sampled from China, Nepal, India, and Thailand, was the 
reciprocally monophyletic sister group of the flavescens-complex. The latter included 
four regional clades: 1) Thailand and Himalaya, 2) China, 3) Malaysia, and 4) the 
Philippines and Taiwan. 

The unique haplotype COI dataset was used to estimate the input tree for the 
species-delimitation analysis. The nucleotide substitution model of this dataset was 
GTR+I (BIC = 3877). Species delimitation using PTP returned each of the nine spe- 
cies as a separate entity, with Bayesian support values between 0.14—1.00 (Fig. 6). In 
contrast, the GMYC lumped the flavescens-complex and B. rotundiceps, as well as the 
B. modestus and B. wangae species pair (Fig. 6). All analyses recovered the flavescens- 
group, including B. flavescens s. str., B. baguionensis, B. bakeri, B. rufoflavus, and B. ta- 
hanensis, as a single species. 


Biogeographic analysis 


The phylogenetic analysis was based on three loci (see supplement for missing data in 16S 
and PEPCK; Suppl. material 6) and 26 terminals, including the outgroup. Tree searches 
were conducted under partitioning and separate model choice for COI (575 bp), 16S 
(551 bp), PEPCK introns (483 bp), and PEPCK exons (369 bp). ‘The resulting maxi- 
mum clade credibility tree showed a deep separation of the B. flavescens/B. rotundiceps 
clade from their closest relatives, B. pratorum and B. ardens, estimated at 6 +/- 2 Ma 
(Fig. 7). Within the flavescens-complex the early split in was occupied by the Malaysian 
group. [he extant members of the flavescens-complex diversified during a time interval 
approximately 1-2 Ma. However, there is some uncertainty about relationships among 
the four regional groups, as in the COI tree the Thai and Himalayan groups were sister 
to the China, Malaysia, the Philippines and Taiwan groups. 

Although DEC model was the best fit model for the BioGEOBEARS analysis in 
this study (the highest AlICc wt; Suppl. material 7), we selected DIVALIKE+J model 
instead (the second high AlCc wt; Suppl. material 7). The DIVALIKE+J model sup- 
ports widespread founder-event speciation at nodes (Yu et al. 2020) which fits the 
distribution changes of bumblebees (Williams et al. 2020, 2022b). 

Results from S-DIVA and DIVALIKE+] showed that there was uncertainty in ances- 
tral areas of B. flavescens and relatives, with more than one possible ancestral area suggested 
in both scenarios (Figs 8, 9). There were only four nodes (III, XI, XIV, and XVI; Fig. 8) 
for which we identified a single ancestral area in S-DIVA, whereas the DIVALIKE+] 
result (Fig. 9), only node XI showed a single ancestral area. The results suggested that 
the ancestral area for the MRCA of the pratorum-group was likely to be in area covered 
by Europe, Central Asia, and the Himalaya (ABF) in S-DIVA, and Central Asia (B) in 
DIVALIKE+] with low probability. This group diversified in the late Miocene. 
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Figure 6. Bayesian inference (BI) phylogenetic tree based on a 575 bp fragment of COI in MrBayes. 


MCMC chains were run for 10 million generations, sampled every 1000" generation. Burn-in fraction is 


10%. The posterior probabilities are shown above the branches. The tip of the tree shows the sample label 


including the sequence length, a taxon name, an identifier code from the database, and its geographic 


origin. The results of species delimitation methods are shown in the right-hand side columns. The black 


and grey bars within the same columns indicate the same species. 


For the flavescens-complex, the S-DIVA biogeographic scenario illustrated that the 
MRCA of B. rotundiceps and B. flavescens (Node IX; Fig. 8) might have occurred in 
the Himalaya (F) around 1.5 Ma during the Pleistocene epoch. Next, the MRCA 
of B. flavescens in the S-DIVA analysis was inferred to be in the FGH areas, includ- 
ing Himalaya, Thailand, and Peninsular Malaysia (Node X; Fig. 8) at around 1 Ma. 
However, the DIVALIKE+J biogeographic scenario suggested that the group might 
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Figure 7. The maximum clade credibility tree of B. flavescens and their closely related species, recon- 
structed with *BEAST from gene trees for the COI, 16S, and PEPCK genes and calibrated using estimate 
time from Hines (2008). The MCMC chains were run for 500 million generations, sampled every 5000" 
generation. Burn-in fraction is 20%. The posterior probabilities are shown under the nodes. Blue bars 


represent the 95% confidence limits on the estimate dates of divergence. Bombus lepidus is the outgroup. 


have originated in Malaysia (H), but with low probability. The DIVALIKE+] showed 
higher support for a Peninsular Malaysian origin (Node X; Fig. 9). 

The diversification of the B. flavescens populations began after 1 Ma, resulting 
in three distinct clades: 1) rufoflavus+tahanensis from Malaysia (Node XI, Figs 8, 
9); 2) flavescens s. str. from Thailand, Himalaya and China (Node XIII, Figs 8, 9); 
3) baguionensis+ bakeri+flavescens s. str. from Taiwan (Node XV, Figs 8, 9). Accord- 
ing to our corridor-dispersal model diagram (Fig. 4), the S-DIVA showed that the 
B. flavescens used all possible corridors between E-K, except the corridor H—K, the 
Sundaland route, whereas the DIVALIKE+J included the corridor H—K as a permit- 
ted route. 


Morphological study of Bombus flavescens ocelli 


The carpenter bees, Xylocopa tranquebarica and X. myops, had a significantly larger 
median ocellus than the bumblebees (Fig. 10). Bombus irisanensis clearly showed the 
smallest median ocellus and B. flavescens, B. rotundiceps, and B. montivagus exhibited 
similar median ocellus sizes (Fig. 10). 
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Figure 8. Biogeographic scenarios for the flavescens-complex and its close relatives by S-DIVA analysis. 
The input tree is reconstructed with *BEAST from gene trees for the COI, 16S, and PEPCK genes and 
calibrated using estimate time from Hines (2008). The internal nodes are represented in Roman numerals. 
The colours of the nodes represent the possible ancestral areas. The area codes used are given in Table 5. 


The most likely ancestral areas are represented above the nodes with the probabilities below. 


For the flavescens-complex (Fig. 11), the Shapiro-Wilk normality test suggested 
that the MOW was not normally distributed (W = 0.93, p-value < 0.05), whereas the 
MOW: ID and MOW: HW did not differ significantly from a normal distribution 
(W = 0.99, p-value = 0.28 and W = 0.99, p-value = 0.45 respectively). There was a 
significant difference within the MOW of the flavescens-complex (Kruskal-Wallis test, 
chi-squared = 36.217, p-value < 0.05; Dunn’s test, p-value < 0.05), including the taxa 
baguionensis-rufoflavus, baguionensis-tahanensis, Himalaya-rufoflavus, and Himalaya- 
tahanensis. For the MOW: ID, nine pairs were significantly different (ANOVA, F-val- 
ue = 6.781, p-value < 0.05; Tukey test, p-value < 0.05): bakeri-China, rufoflavus-Chi- 
na, tahanensis-China, bakeri-Himalaya , rufoflavus-Himalaya, tahanensis-Himalaya, 
Taiwan-bakeri, Taiwan-rufoflavus, and Taiwan-tahanensis. In addition, only two pairs 


The oriental bumblebee Bombus flavescens 523 


@ 4- Europe @ 6-Thaitana Pie (J) baguionensis 
0.50 
@ 8-casia @ b- Malaysia kal (I) flavescens (Taiwan) 
© c-NAsia+NChina @ 1- Taiwan oas a XV 
@ D - Korea + Japan @ J- Luzon (K) bakeri 


7 4 | 0.30 
= : as @ «Negros (E) flavescens (China) 
- Himalaya 


= @ Xill , 
049 (F) flavescens (Himalaya) 


L@xiv 
0.60 a5 (G) flavescens (Thailand) 
P (H) rufoflavus 

H 


0.04 


Xl 
at (H) tahanensis 
0.23 vill (EFG) rotundiceps 


(A) pratorum 
A amy "7 (D) ardens 


(CD) modestus 


0.53 vl 
B Ad . (E) wangae 


(A) pyrenaeus 
7 Ty (B) biroi 
0.03 
(F) nursei 
(EF) /epidus 

8 7 6 5 4 3 2 1 O Millions of years ago 
Miocene Pliocene | Pleistocene 

Figure 9. Biogeographic scenarios for the flavescens-complex and its close relatives by BioGeoBEARS with 
DIVALIKE+J analysis. The input tree is reconstructed with *BEAST from gene trees for the COI, 16S, and 
PEPCK genes and calibrated using estimate time from Hines (2008). The internal nodes are represented in 
Roman numerals. The colours of the nodes represent the possible ancestral areas. The area codes used are 


given in Table 5. The most likely ancestral areas are represented above the nodes with the probabilities below. 
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Figure 10. Box plots showing the median ocellus width (MOW), the ratio between MOW and inter- 
tegular distance (MOW:ID), the ratio between MOW and the head width (MOW:HW) of B. irisanensis, 
B. rotundiceps, B. flavescens, B. montivagus taxon maxwelli, Xylocopa tranquebarica, and X. myops. 
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Figure I |. The box plots of the median ocellus width (MOW), the ratio between MOW and intertegular 
distance (MOW-ID) and the ratio between MOW and head width (MOW:HW) among Bombus flavescens 
population: Thailand, China, Himalaya, Taiwan, baguionensis, bakeri, rufoflavus, tahanensis. 


of taxa were significantly different in MOW:HW (ANOVA, F-value = 3.701, p-val- 
ue < 0.05; Tukey test, p-value < 0.05), rufoflavus-Himalaya, and tahanensis-Himalaya. 


Discussion 


Species status of the flavescens-complex 


The status of species within the flavescens-complex taxa has been debated for nearly a 
century. Theodore Frison (1895-1945), an American entomologist wrote in his work 
in 1934, “This species of bumblebee [B. flavescens| has been a problem to most persons 
who have attempted to determine...” (Frison 1934). Our results for the COI-barcode 
tree (Fig. 6) support that the flavescens-complex is a monophyletic group, including 
the populations in China, the Himalaya, and Southeast Asia. The PTP and GMYC 
results group the flavescens-complex together as one species. This demonstrates that 
B. flavescens s. l. includes high intraspecific variation in colour pattern across its distri- 
butional range. The species-delimitation analyses confirm the conspecific status of the 
various geographical and morphological variants named as B. baguionensis, B. bakeri, 
B. tahanensis and B. rufoflavus by previous authors. Although this study did not in- 
clude the taxon mearnsi from Mindanao, the taxon mearnsi also is expected to be a 
synonym of B. flavescens as morphologically suggested by Pittioni (1949), similar to the 
taxon bakeri from Negros, another population from the southern Philippines islands. 
There is a discrepancy between the number of species recognised by PTP and 
GMYC. The GMYC analysis lumped 1) B. flavescens and B. rotundiceps 2) B. modestus 
and B. wangae, as single species, whereas PTP split them into four species. GMYC 
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requires an ultrametric tree, which distorts the tree and branch lengths (Fujisawa and 
Barraclough 2013). For this, the closely related sister taxa which establish their own 
short clade might be identified as the same species. In contrast, the PTP keeps the orig- 
inal tree shape, and should be considered more reliable between the two approaches 
(Williams 2021), in this instance. 

Nevertheless, when we investigate the morphological evidence of these four taxa, 
their morphological characters are unique (Williams 1998; Williams et al. 2009, 2010). 
The males of B. wangae and B. modestus are distinct, especially in their genitalia, which 
the gonostylus of B. modestus is broader than B. wangae. For the B. flavescens and B. ro- 
tundiceps, their male genitalia are also different: the gonostylus of B. flavescens is broad 
with round inner margin, whereas the gonostylus of B. rotundiceps is narrower with 
straight inner margin. In this case, the PTP results show congruence with the morpho- 
logical characters whereas GMYC does not. Then, we suggest that the PTP analysis cor- 
roborates more consistently morphologically recognisable species for these bumblebees. 
For further study, if fresh male specimens of B. flavescens are available, chemical evidence, 
including, CLGS, is recommended as an alternative approach for defining the species 
recognition, whether it shows agreement with the molecular delimitation in this study. 

Fresh specimens of B. flavescens and their relatives are not available from enough 
samples because of the rarity of these bees. Bumblebees in subgenus Pyrobombus are 
active particularly early in the year (Williams et al. 2014), so subgenus Pyrobombus 
sampling can only be conducted in a relatively short period. The accessibility of their 
habitats in Southeast Asia is also relatively limited and no records have been reported 
during recent decades, especially for B. flavescens from Gunung Tahan, Malaysia, and 
from Negros and Mindanao in the Philippines. Moreover, although the flavescens-com- 
plex is widely distributed, recorded in at least eleven countries, it is difficult to obtain 
specimens from some countries, because of limits imposed by local or international 
regulations. Therefore, old museum specimens are the best available option for clarify- 
ing genetic relationships within this group. However, these specimens were collected 
between 1980 and late 1990. The DNA of those specimens has become degraded 
through time (Paabo 1989). It is difficult to extract the DNA from old specimens 
with standard protocols. However, numerous DNA extraction techniques for old in- 
sect specimens have been developed (Gilbert et al. 2007; Thomsen et al. 2009; Ballare 
et al. 2019). Next generation sequencing would help to obtain DNA sequences from 
historical specimens (Nazari et al. 2016; Prosser et al. 2016; Sproul and Maddison 
2017; Call et al. 2021), but a high DNA quantity is required (Goodwin et al. 2016). 
According to our results, the DNA concentration of the samples in this study was not 
high enough for the sequencing requirement. The newly designed primers in this study 
were successful for the flavescens-complex and B. rotundiceps. The primers can amplify 
the DNA from specimens up to 70 years old. From our results, fresh specimens or 
pinned specimens less than 20 years old are recommended for use as DNA sources as 
the first option. Museum specimens are an alternative way if fresh specimens are not 
available. However, using museum specimens requires additional molecular processes, 
for example, extra-sterilisation of equipment and laboratory space, and conducting 
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negative extraction and PCR controls. ‘This is the first time that at least the partial COI 
barcodes of the flavescens-complex populations from the Philippines from both North 
and South islands have been retrieved. This information is invaluable because this is 
the key to resolving cryptic status among this group, but it also provides a workflow 
for approaching other difficult bee groups for which only old specimens are available. 


Orange pattern bumblebees in Peninsular Malaysia 


In the Cameron Highlands, Malaysia, several orange pattern bumblebees have been 
recorded, including the rufoflavus form of B. flavescens (Fig. 2J) and the maxwelli form 
of B. montivagus. Their fulvous hair and sclerite colour might relate to the nocturnal 
lifestyle of hymenopterans, similar to the nocturnal carpenter bees in subgenus Nyc- 
tomelitta (Williams 2007). In this study, the ocelli of B. flavescens populations were 
similar in size (Fig. 11) which did not match the hypothesis of enlarged ocelli, whereas 
the nocturnal carpenter bees have much more distinctively enlarged ocelli (Fig. 10). 
However, in the context of the variation among B. flavescens populations, the median 
ocellus of both Malaysian populations, B. flavescens taxon rufoflavus and taxon tahan- 
ensis, were significantly larger than the populations at higher latitudes, including from 
Himalaya and Luzon Island (Fig. 11). The results showed a gradient of ocellar size 
within the B. flavescens populations from high latitude (smaller ocellus) to low latitude 
(larger ocellus): MOW, MOWID, and MOWHW differed significantly between low 
latitude (0°—20°N) and high latitude (>20°N) group (MOW: Mann-Whitney U test, 
p-value = 0.037; MOWID: ANOVA, F-value = 30.49, p-value < 0.05; MOWHW: 
ANOVA, F-value = 8.573, p-value < 0.05). 

Although there is a significant difference between low-latitude and high-latitude 
populations of B. flavescens, intraspecific variation of ocellar size in bumblebee popu- 
lations has been reported before for B. terrestris (Kapustjanskij et al. 2007), without 
any evidence for nocturnal behavior in this well-known species. ‘This also suggests that 
ocellar size might not be a good diagnostic character to distinguish species in bumble- 
bees; it may be that it is relatively plastic to local selection regimes. ‘The ocelli are visual 
organs that detect polarized light in low-light conditions (Wellington 1974; Roubik 
1989). Ocellar size shows a negative correlation with light intensity (Kerfoot 1967; 
Warrant et al. 2006; Kapustjanskij et al. 2007), and there are a high number of noctur- 
nal and crepuscular bees in tropical or low latitude areas (Dorey et al. 2020). The trend 
seen here towards larger ocellar size in the tropics might be evidence that bumblebees 
are more active in low-light conditions, for example, in the early morning or late after- 
noon, to avoid heat stress, as is suggested for B. breviceps Smith and B. haemorrhoidalis 
Smith (Williams 1991; Thanoosing 2022). Larger ocelli can also be observed in other 
tropical forest bees, for example, stingless bees, because sunlight is filtered by the tree 
canopy, thereby reduced in the understory (Streinzer et al. 2016). 

Numerous light traps were run at night on the Cameron Highlands recently by re- 
searchers (Musthafa et al. 2021) and by tourist services (e.g., Cameron Service: https:// 
cameronservice.blogspot.com/), but no bumblebees have been observed at traps. Con- 
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sequently, a nocturnal or crepuscular lifestyle of the orange pattern bumblebees in the 
Cameron Highlands remains unsupported. 

Apart from orange B. flavescens, another bumblebee, B. montivagus taxon max- 
welli, and a variation of the diurnal hornet Vespa velutina Lepeletier taxon divergens, 
in the same locality, are also orange in colour (Hines et al. 2012; Perrard et al. 2014). 
In addition, not only the orange pattern but also the black hair with a red tail pattern 
of B. flavescens in South China (Fig. 2D) is recognised in the colour variation of V. ve- 
lutina taxon nigrithorax (Perrard et al. 2014). For this reason, the colour patterns of 
B. flavescens might potentially be a result of Miillerian mimicry. 


Biogeography of Bombus flavescens 


Pyrobombus is a bumblebee subgenus in which five species groups are found in the 
Old and New World (Williams 1998). Bombus flavescens and its relatives are in the 
‘pratorum-group. This group is the sister group to the ‘/epidus-group’ (Cameron et 
al. 2007). The groups originated in the Palearctic and Oriental regions, in the late 
Miocene (~11 Ma), and the /epidus-group is entirely restricted to the Oriental region 
(Williams 1998). In this study, biogeographic analyses show that the origin of the pra- 
torum-group is likely to be in the Oriental region, specifically around the Qinghai-Ti- 
betan Plateau (QTP) or Himalaya. Then, during the late Miocene and early Pliocene, 
the diversification of this group occurred. This is coincident with the period when the 
monsoon climate in this area intensified (Zhou et al. 2018), which also accelerated the 
diversification of bumblebee food plants in the Oriental region (Yu et al. 2015; Ma 
et al. 2016; Matuszak et al. 2016). Our results show that there were dispersal events 
between the Oriental (E, FE, G, H, I, J, and K; Fig. 8) and the Palearctic region (A, B, 
C, and D; Fig. 8), for example, the clade of biroi+nursei, pyrenaeus+ modestus+ wangae, 
and pratorum+ardens. \n addition, there are the lineages of B. nursei and B. wangae 
that reinvaded the Oriental Region: B. nursei from the west and B. wangae from the 
east. Bombus nursei prefers subalpine meadow habitats similar to B. biroi (Williams 
1991; Williams 2022). This dispersal pattern between the Oriental and Palearctic re- 
gions during the Miocene-Pliocene can be observed in other bumblebee subgenera, 
for example, in subgenera Mendacibombus (Williams et al. 2018) and Melanobombus 
(Williams et al. 2020). 

Only the group rotundiceps+flavescens has been restricted to the Oriental until now. 
The divergence between B. rotundiceps and B. flavescens was in the Pleistocene (ca 1.5 Ma). 
The climate in Southeast Asia during the late Pliocene and Pleistocene was cooler than 
the present day (Morley 2012). This might have given the opportunity for temperate 
and montane flora to colonise these areas. At the same time, the MRCA of B. flavescens 
dispersed southward into Southeast Asia or Sundaland, following suitable habitat (A; 
Fig. 12). Bombus flavescens is likely to have dispersed to and colonised the Philippines 
islands when the islands or the Pleistocene Aggregate Island Complex (PAIC), including 
Luzon, Negros-Panay or Visayan, and Mindanao, were connected in the LGM (Brown 


et al. 2013). This study supports that B. flavescens dispersed to the Philippines by the 
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A 15Ma Bima 


C 0.5Ma D present 


Figure |2. Biogeographic scenarios of Bombus flavescens through time A the most recent ancestor of B. 
flavescens dispersed through Southeast Asia or Sundaland around 1 Ma (orange) B the Malaysian popu- 
lation was isolated on the highland of the peninsula due to the temperature rising around 1 Ma (blue), 
at that time, the populations in Taiwan and Philippines were connected to the mainland populations 
(orange) C the bridge between the mainland and Taiwan was submerged, the Taiwan-Philippines popula- 
tions were isolated around 0.5 Ma (purple) D at present day, B. flavescens populations distribute in the 
subtropical of Southeast Asia, China and Himalaya (yellow), on the Highlands of Malay peninsula (blue), 
and the islands of Taiwan and the Philippines (purple). 


Taiwan-Luzon route in the same way as B. irisanensis, the other species in the Philip- 
pines. The oscillation of sea level and raising of global temperature in the late Pleistocene 
facilitated the population isolation and higher latitude immigration of B. flavescens (B, 
C; Fig. 12). First, the populations from the North (Luzon) and the South (Negros and 
Mindanao) were separated, because the South population (taxon bakeri) is the sister of 
the North population (taxon baguionensis) and the Taiwan group (Fig. 7). Then, the 
populations on the North of Philippines and Taiwan were isolated when the sea sub- 
merged the bridges, so that the populations in Southeast Asia were captured in montane 
refuges as biological islands until the present, similar to the populations on the Cameron 


Highlands and Gunung Tahan mountain in Malaysia (D; Fig. 12). 


Biogeography of bumblebees in Southeast Asia 


Many lineages of bumblebees in Southeast Asia might are hypothesised to have origi- 
nated around the QTP, then dispersing into the region via the Himalaya-Hengduan 
corridor (Williams 1985; Hines 2008; Williams et al. 2022). The bumblebee fauna 
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at the Northern border and highlands of the Southeast Asian region is similar to the 
fauna of the Himalaya and Southern China, adjacent neighbours both spatially and 
genetically (Williams et al. 2009; Williams et al. 2010; Hines and Williams 2012; 
Streinzer et al. 2019; Williams 2022). The subgenera Orientalibombus Richards and 
Alpigenobombus Skorikov are mainly restricted to the North of the Southeast Asian 
region, whereas the subgenera Megabombus Dalla Torre, Melanobombus Dalla Torre, 
and Pyrobombus are more widely distributed throughout the region (Williams 1998). 

The long-faced subgenus Megabombus crown group might have originated around 
13 Ma (Hines 2008). Then, two independent lineages of subgenus Megabombus di- 
verged in Southeast Asia approximately between 4.25—1.5 Ma, based on molecular dat- 
ing (Huang et al. 2015a). The lineages are 1) a trifasciatus-group: B. montivagus, B. albo- 
pleuralis Friese, B. burmensis (Skorikov) in the North of Southeast Asia and the Malay 
Peninsula, and 2) a senex-group: B. irisanensis on Luzon Island, in the Philippines. 
However, due to a lack of genetic information on B. senex Vollenhoven and B. melan- 
opoda Cockerell, the phylogenetic relationship between mainland and Sumatran subge- 
nus Megabombus remains relatively unresolved (Huang et al. 2015a). Consequently, the 
biogeography of Sumatran subgenus Megabombus bumblebees remains unclear. 

Nonetheless, the link between the Southeast Asian mainland and the Indonesian 
islands can be explained by the phylogenetic relationship of two subgenus Melanobom- 
bus sister species, B. eximius Smith and B. rufipes Lepeletier of the rufipes-group. The 
MRCA of the rufipes-group lineage diverged from the other subgenus Melanobombus 
lineages around 16 Ma and was distributed in the Himalaya and QTP (Williams et al. 
2020). Speciation within the rufipes-group is likely to have occurred around 1 Ma in 
the Pleistocene epoch, after the Sundaland land bridge was submerged, isolating part 
of the rufipes-group on Sumatra and Java (Williams et al. 2020). 

Southeast Asian bumblebees in the subgenus Pyrobombus are also distributed both 
on the mainland and on the islands of the Philippines, for example, B. flavescens. There 
is no record of subgenus Pyrobombus on the islands of Indonesia and adjacent areas. A 
member of subgenus Pyrobombus had been recorded on the Andaman Islands, named 
B. andamanus Gribodo (Gribodo 1882). However, this is a mislabeled specimen of a 
Nearctic bumblebee, B. bifarius Cresson (Williams 1998). 


Conclusion 


Genetic information proved crucial for the study of bumblebees in Southeast Asia. 
This is the first gene-based study to address the taxonomic status of B. flavescens, which 
is highly variable in hair colour. This study confirms that the populations of B. flave- 
scens on the Asian mainland and on the islands are parts of the same species. Despite 
this, a trend among B. flavescens populations can be observed towards larger ocelli at 
lower latitudes. This might only reflect local selection based on foraging in the low- 
light conditions within dense forest, or in more early morning or twilight activity in 
warmer environments. Bombus flavescens originated in the Himalaya and dispersed to 
Southeast Asia during the Pleistocene. Its constituents, various regional colour forms, 
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diversified through an allopatric divergence process. Bombus flavescens is a useful model 
for studying the biogeography of bumblebees in Southeast Asia, many of which are less 
known. Nevertheless, more genetic information is required to investigate the conserva- 
tion of endemic populations of B. flavescens. 
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